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Abstract 

We study the high-velocity regime mode-I fracture instability when small microbranches start 
to appear near the main crack, using large scale simulations. Some of the features of those mi¬ 
crobranches have been reproduced qualitatively in smaller scale studies (using 0(1O 4 ) atoms) on 
both a model of an amorphous materials (via the continuous random network model) and using 
perturbed lattice models. In this study, larger scale simulations (0(1O 6 ) atoms) were performed 
using multi-threading computing on a GPU device, in order to achieve more physically realistic 
results. First, we find that the microbranching pattern appears to be converging with the lattice 
width. Second, the simulations reproduce the growth of the size of a microbranch as a function of 
the crack velocity, as well as the increase of the amplitude of the derivative of the electrical resis¬ 
tance RMS with respect to the time as a function of the crack velocity. In addition, the simulations 
yield the correct branching angle of the microbranches, and the power law governing the shape of 
the microbranches seems to be lower than one, so that the side cracks turn over in the direction of 
propagation of the main crack as seen in experiment. 
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I. INTRODUCTION 


The study of the physics of brittle fracture has been very fruitful in the last two decades ID- 
13]. New experiments have shown various new features of dynamic fracture, focused in Mode- 
I (tensile) fracture using amorphous brittle materials [3H7]. In particular, the experiments 
have shown a sharp transition between the regime of steady-state cracks and the regime of 
unstable cracks mn- In steady state, where the driving displacement, A (of order Aq, 
the Griffith criterion), is sufficiently small, a single crack propagates in the midline of the 
sample, reaching a steady-state velocity (which is of order the Rayleigh surface wave speed, 
cr). Increasing A results in an increased steady-state velocity, yielding a v(A) curve, until 
a specific critical point. Increasing the diving displacement further, beyond this critical 
point, we enter the unstable regime, where small microbranches start to appear nearby the 
main crack. The experiments have shown that above the critical point, the size of the aver¬ 
age microbranch (which is log-normal distributed) increases rapidly with the crack velocity, 
measured via the slope of the electrical resistance of a conductive layer that is pasted on 
the sample. The electrical resistance slope exhibits oscillations whose amplitude increases 
rapidly as well. Increasing the driving displacement furthermore, the small microbranches 
become large microbranches, creating a complex fracture pattern, and finally, creates mac¬ 
robranches mm. 

Some of the new experimental findings could not be explained via the classical theoretical 
approach for fracture mechanics, the linear elasticity fracture mechanics (LEFM) theory [9]. 
For example, several predictions of LEFM for the critical velocity, such as the studies of 
Yoffe [10] or Eshelby ED predicted a single universal critical velocity much higher than 
that seen in some of the experimental studies, such as in PMMA pQ. Also, experiments 
have found various material-dependent features such as the terminal velocity which the 
crack manages to propagate as well as the critical velocity for macro-branches [I]. As far 
as the micro-branching critical velocity, the question of universality is debatable H2. but 
anyhow, the velocities are much smaller that the theoretical predictions. For a review, see 
the Introduction in [13]. However, the basic reason for the failures of LEFM is that the 
basic equations of LEFM yield a singularity of the stresses near the crack tip [9], and thus, 
the zone nearest to the crack’s tip (the process zone) cannot be modeled via LEFM. 

The failure of LEFM, caused by this singularity, gave rise to an interest in discrete lattice 
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models and simulations [HH23]. In this kind of model, the basic length-scale which is the 
lattice scale, prevents the singularities that appear in the continuum approach. The lattice 
models were successful in reproducing the behavior of steady-state cracks puma EH, includ¬ 
ing a material dependency of the v(A) curves, in that it depends on the specific parameters 
of the inter-atomic potential mm- Moreover, the lattice models predict a certain criti¬ 
cal point beyond which the steady-state solution becomes linearly unstable [ISM]. In the 
simulations, this is exactly the point which the crack stops propagating along the midline of 
the sample and some additional bonds, not along the midline, start to break [TUI IS1 EDI f2T| • 
However, especially in mode-I pure lattice simulations, the post-instability behavior of the 
lattice models do not match the experiments, neither qualitatively or quantitatively [201 El]. 

One recent approach to overcome these difficulties has been to turn to a more realistic 
model for an amorphous material, the continuous random network (CRN) model [25] . The 
continuous random network was suggested first by Zachariasen for describing amorphous 
material [2U], and specific effective algorithms (using Monte-Carlo techniques) for generating 
the CRN were offered in [251 E3 EE]. Recent accurate 2D experiments using transmission 
electron microscopy (TEM) in 2D silica on the structure of this amorphous material were 
reproduced to excellent accuracy using the Zachariasen model [29, [30]. After generating 
an amorphous CRN sample, molecular dynamics simulations were found to yield many 
of the important qualitative features of mode-I fracture experiments in amorphous brittle 
materials [25]. The simulations showed the birth of microbranches growing nearby the main 
crack, the increase of the size of the microbranches as a function of the driving displacement 
(or of course, the crack velocity), and the growth of the amplitude of the derivative of the 
electrical resistance with respect to the time as a function of the crack velocity. 

Another direction recently examined was that of perturbed lattice models [13], which 
exhibited behavior similar to that of the CRN model, including the main features mentioned 
above. However, these simulations (both the CRN and the perturbed lattice simulations) 
suffered from a significant level of numerical noise, since each microbranch contained only 
a few tens of broken bonds at most. Thus, the statistics that concerns the most interesting 
physics, that of the branches, was quite poor. 

These intriguing results, performed on limited size systems (0(5 ■ 10 4 ) atoms) provide 
strong motivation to conduct larger simulations both to reduce the overall noise level and to 
get closer to at least a mesoscopic system where scaling behavior might set in. The goal was 
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to achieve at least a two magnitude increase in size, i.e., simulations of the order of millions 
of atoms (0(5 ■ 10 6 ) atoms), which necessitated using parallel (multi-threading) computing. 
In this work we study mode-I fracture via large scale simulations, with a particular focus on 
the properties of the microbranches that could not be studied in the previous, limited size, 
studies. In addition, recently experiments have been performed with gels, whose three orders 
of magnitude slower sound speed (compared to the classic brittle materials like PMMA or 
glass) enables direct visualization by means of moderately fast video cameras, resulting in 
clear snapshots of the crack tip Ha ca mu. Using the larger scale simulations we can now 
compare the crack’s tip shape, both on steady-state cracks in a strip geometry, and especially, 
near the origin of the microbranching instability. 

The models are presented in Sec. [TT[ In Sec. |III| we perform some basic checks of our 
models, confirming that the transverse size of the microbranch zone ( 5 v/w) decreases with 
increasing sample width W, for a given scaled driving displacement A/Aq (in the exper¬ 
iments, using macroscopic sample sizes, the microbranching region width is much smaller 
than the sample width and the dynamics of the fracture is not effected by the sample edges). 


Next, in Sec. IV we present the quantitative results concerning the birth and the growth of 
the microbranches, and their physical features. A short discussion is presented in Sec. |V) 


II. MODEL AND MAIN METHODOLOGY 

The simulations presented in this study are divided generally into two kinds. The first 
one uses the continuous random network (CRN) model (to model an amorphous material), 
and the second employs a perturbed (honeycomb or hexagonal) lattice model. Both models 
were described in depth in [25] and ra respectively. We will review them here shortly. 

The CRN model reproduces various structural features of real brittle amorphous ma¬ 
terials (though 2D, in this study) such as amorphous silicon or silica [29] [30], 32], and a 
3D-extension of this model should reproduce the real behavior of fracture. The perturbed 
honeycomb lattice is the ordered phase of the CRN, as discussed at length in [13] and shares 
similar features and results (though with less noise). In addition, in [29l 130] there is clear 
experimental 2D evidence that 2D ordered materials share the structural features of the 
perturbed honeycomb lattice. Also, the perturbed hexagonal lattice is a generalization of 
the hexagonal lattice that was used in many studies to study dynamic fracture (for exarn- 
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pie [H 10, [23]), and facilitates comparison with this body of work. 


1. Generating the Continuous Random Network for modeling amorphous material 


We generated two-dimensional CRN’s by a 2D-analogue [ 25] of the WWW algorithm [2?, 
[28]. The potential that was used in the construction of the CRN included both a 2-body 
central force and a 3-body bond-bending force: [33, 33]: 


Etot - Y 
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where \rfj\ is the radial distance between each pair of nearest-neighbor atoms and a t j = 
ao = 4 is a constant lattice scale (in contrast to the perturbed lattice model). k r and kg are 
the radial and the azimuthal (3-body) spring constants, respectively, cos I s the cosine 
of the angles between each set of 3 neighboring atoms, defined of course by: 


COS 9i : j : k 


^i,j * k 


( 2 ) 


where i is the central atom and (j, k) are its two neighbors. Q c = 2?r /3 (characterizing 
a honeycomb lattice). We start from a pure honeycomb lattice, randomize large number 
of bonds and perform a Monte-Carlo procedure, wherein each cycle we switch two bonds, 
calculating the optimal positions of the atoms in the near zone of the switched bonds to 
determine the change of energy so as to decide whether to accept the switch. Finally, we 
get a CRN that looks like Zachariasen’s patterns |25] 26] and the TEM snapshots of the 2D 
amorphous Silica 2TJ. SI]- For a in-depth discussion, see [23j. 


2. Generating the perturbed lattice 

Here we start with a perfect honeycomb lattice and randomize the lattice scale of each 
“bond”, a t ,f 

c^i.j (l T 0) i 1)2,... ,n a toms)J £ <A/"(i) (3) 

where e l:J G [—b,b], and b is constant for a given lattice, and in this work ranges between 
0^6^ 0.1, a 0 = 4. A f[i) refers to the nearest-neighbors of site i. For a detailed discussion, 
see [ 13] . 
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3. The equations of motion 


In our model, between each two atoms there is a piece-wise linear radial force (2-body 
force law) of the form: 

flj = krk'ij (| f id | - a id )fj, h (4) 

where: 

K,j = 0 H (e- \r itj \). (5) 


The Heaviside step function 6 # guarantees that the force drops immediately to zero when the 
distance between two atoms |r)j| reaches a certain value e > (the breaking of a “bond”). 
In this work we set £ = do + 1. We describe here brittle materials with an extreme sharp 
transition from linear response to failure, though, in reality the failure is always somewhat 
smoother. Thus, alternatively to Eq. [5j we can use a smoother force law, which instead of a 
sharp failure at |f)j| = £, has a more realistic smooth transition wherein the force law drops 
to zero of the form [18, 2D1 l2Tj : 


. 1 + tanh[a pot (e - fjj)] 


( 6 ) 


1 + tanh(a pot ) 

where a pot is the smoothness parameter, such that when cc pot —>• oo the force law reverts to 
the piecewise linear force law. The effect of a pot on the fracture feature was investigated 
previously da edi ei], and is reproduced in this paper. The results in this paper refer to 
the piecewise linear model, unless mentioned otherwise. 

In addition there is a 3-body force law that depends on the cosine of each of the angles, 
acts on the central atom (atom i) of each angle, and may be expressed as: 
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while the force that is applied on the other two atoms (atoms j, k) may expressed as: 
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Of course, the forces satisfy the relation: ffu k \ = ~{fju k ) + f k (ij ))■ The 3-body force law 
drops immediately to zero when using a piecewise linear force law when the bond breaks 
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(Eq. g, or may be taken to vanish smoothly, using Eq. [6j In a honeycomb lattice there are 
three angles associated with each atom and in the hexagonal lattice there are six of them 
(we note that in the hexagonal lattice this choice is a little bit arbitrary since there are in 
general additional optional angles for each atom). There is a certain preferred angle 6c for 
which the 3-body force law vanishes (in the honeycomb lattice we set 6 C = 2 Vs and in the 
hexagonal lattice we set 6c = V 3 )- 

In addition, it is convenient to add a small Kelvin-type viscoelastic force proportional to 
the relative velocity between the two atoms of the bond [171120] : 

9i,j = v(vij-f i j)k' i / itj} (9) 

with 7] the viscosity parameter. The viscous force also vanishes after the bond is broken, 
governed by k\y The imposition of a small amount of such a viscosity acts to stabilize the 
system and is especially useful in the relatively narrow systems simulated herein. 

The set of equations of motion of each atom is then: 

m-ai= (./J, + 9i,j) + ^2 + ^2 fjfak)’ ( 10 ) 

j£3p nn j,k£3pnn j£6p nn 

In this work the units are chosen so that the radial spring constant k r and the atoms mass 
rrii is unity. We note that numerical measurement of the Young’s modulus of the hexagonal 
lattices for kg = 0 yields the well-known analytical expression 2/\/3 k r and is twice as big 
(E ps 2 k r ) with kg/k r = 10 and the Poisson’s ratio is is = 1/3 as reported in many previous 
works [35j [36]. In the honeycomb lattice and the CRN (with k e /k r = 1), the Young’s 
modulus is approximately E ~ 0.86fc r while the Poisson’s ratio remains similar. However, 
the value of k r is not crucial (in this study), since the results in this paper are in normalized 
units (A/A g , v/c r ). 

After relaxing the initial lattice, we strain the lattice under a mode-I tensile loading with 
a given constant strain corresponding a given driving displacement ±A of the edges (the 
top and bottom rows are held tight and do not allow transverse displacement) and seed the 
system with an initial crack (The left boundary condition is also held in a pure “cracked” 
state). The crack then propagates via the same molecular dynamics Euler scheme using 
Eqs. lUfToj We note that the calculation of A G is set under the equality of the energy in 
the uncracked uniform strain, to the energy needed for breaking the midline bonds in the 
sample (for example, see [T7J for the case of square lattice). We parameterize our results in 


7 


this paper using the normalized quantity A/Ac, but of course A/A q = Kj/Kic (the stress 
intensity factor normalized to the Griffith value) [Ij. 

4- Parallelization by GPU computing 

As mentioned in the Introduction, the major innovation of this work, compared to our 
previous studies, is the use of large scale simulations. The previous studies of the amorphous 
(CRN) model [25] and the perturbed lattice model pjjj used approximately 50, 000 particles, 
in this study we wished to use approximately 5, 000, 000 particles. These kind of simulations 
cannot reasonably be performed by a singe CPU, and thus force us to use multi-thread 
computing. We choose to use GPU computing, parallelizing the code via CUDA [STUM]. 
This kind of programing forces the programmer to use the different levels of memory care¬ 
fully [38], which makes possible achieving an acceleration up to ~ 100 faster than a regular 
C code using a single CPU. Beside the benefit of getting the results in a given system much 
faster, the main benefit is the possibility to run large scale simulations, which was prohibitive 
before. This tool makes possible the simulation of millions of atoms in reasonable simulation 
times. 

Our model consists of several modules, each one of which needs to be re-written in CUDA. 
Both amorphous and lattice models use a molecular-dynamics module for the fracture sim¬ 
ulations that must be parallelized. In addition, for the CRN, the Metropolis Monte-Carlo 
algorithm for generating the CRN needs to be parallelized. Furthermore, the electrical resis¬ 
tance simulations, which we use to determine the crack velocity [25], is solved by a nonlinear 
Laplace solver that needs to be parallelized too. The electrical resistance is calculated via 
the method that was used in ESI, by solving the nonlinear Laplace on a grid of resistors 
(each broken bond, in the main crack and in the microbranches is taken into account) [39] . 
In Appendix [A] we discuss about the implementation of the different modules, the param¬ 
eters that were used and the degree of acceleration achieved for the different modules with 
the GPU using CUDA. 

The various sized lattices we use contain: 

• 162 ■ 272 « 45, 000 (N = 80 in the Slepyan model notation) atoms for the honeycomb 
lattice and 162 • 408 ~ 65, 000 atoms for the hexagonal lattice, which we call / = 1 
(Factor=l). This size is equal to that used in our previous studies, [25] and [13] . 


• 486-816 ~ 400, 000 (N = 240 in the Slepyan model notation) atoms for the honeycomb 
lattice and 486 • 1224 « 600, 000 atoms for the hexagonal lattice, which we call / = 3 
(Factor=3). 

• 1458 • 2448 ~ 3, 600, 000 (N = 720 in the Slepyan model notation) atoms for the 
honeycomb lattice and 1296 • 3264 « 4,200,000 (N = 640 in the Slepyan model 
notation) atoms for the hexagonal lattice, which we call / = 8 — 9 (Factor=9 in the 
honeycomb and CRN lattices and Factor=8 in the hexagonal lattice). 

In Appendix [B] we present a brief discussion regarding the results for the CRN using the 
parallel GPU Monte-Carlo algorithm. The GPU algorithm reproduced the results of the 
CPU code, in particular the agreement ESI of the radial distribution function with that 
experimentally determined [32] for amorphous silicon. 


III. OVERVIEW OF THE SIMULATION RESULTS 


In Figs. S a ) to |3 c) we present the fracture pattern of the broken bonds for the small 
size perturbed honeycomb lattice (that was used before in [T3J, called here / = 1), for 
the intermediate size lattice (/ = 3), and for the large lattice (/ = 9), for three values 
of the driving displacement: small (A/Ac = 2.8), intermediate (A/Ac = 3.4), and large 
(A/Aq = 4). The fracture patterns are plotted in the x — y plane and depict the full 
simulated sample, when x and y have the units of do, the lattice scale. In the large driving 
simulation for / = 1 reported upon in our previous work, the microbranches reached the edge 
of the sample. In addition, in Fig. [3(a) we present fracture patterns using the amorphous 
CRN model (that was used before in [25]) with A/Aq = 3.6 for the different sizes of lattices, 
and in Fig. [3(b), fracture patterns using a perturbed hexagonal lattice. 

In the Figures we can immediately see the benefit of the larger scale simulations; the 
noisy fracture patterns that were obtained using / = 1 (the upper pattern in each figure), 
transform to the smoother and more physical-like patterns at / = 8 — 9. In Fig. [3(b) we can 
see clearly the curved power-law shape of the microbranches (for a quantitative discussion, 


see Sec. IV C) 


A closer look reveals an important point: For a given driving displacement, the relative 
width of the fracture pattern decreases with the increase of the lattice size. This is crucial 
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Honeycomb A/A G =2.8 Honeycomb A/A G =3.4 




(a) X (b) 

Honeycomb A/A Q =4 



FIG. 1. (Color online) (a) The microbranching pattern in a perturbed honeycomb lattice for 
A/A q = 2.8 and 7] = 2 for / = 1 in the upper curve, for / = 3 in the intermediate curve and for 
/ = 9 in the lower curve, x and y have the units of do, the lattice scale, and each figure depicts the 
whole sample (except for the the initial seed crack that extends from starts from the left edge of 
each crack pattern to x = y = 0). (b) The same for A/ Aq = 3.4. (c) The same for A /Ac = 4.0. In 
(b) we define 5y, the width of the microbranching region, as the difference between the maximum 
and minimum y’s of broken bonds. 
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FIG. 2. (Color online) (a) The microbranching pattern using the continuous random network 
(CRN) model for A/Ac = 3.6 and y = 2 for / = 1 in the upper curve, for / = 3 in the 
intermediate curve and for / = 9 in the lower curve for continuous random network (CRN) model 
(b) The same for perturbed hexagonal lattice for A/A q = 1.7 and y = 0.25. 


since otherwise the branching pattern in a macroscopic material would be macroscopic as 
well, against the evidence of the experiments [5]. In Fig. [3j we can see the relative fracture 
zone width s ii/w, when 5y is the width of the microbranching pattern (see Fig. 0b) for a 
pictorial explanation), defined as the difference between the maximum and minimum y 's of 
broken bonds, and W is the sample width. For any given value of a /a g the normalized width 
of the microbranching pattern decreases with the lattice width. This effect in seen clearly 
in the perturbed honeycomb lattice and in the CRN lattice, and in a more moderate way 
in the hexagonal perturbed lattice. This result is crucial; if the lattice models are physical, 
then when increasing the lattice size, the relative fracture zone must decrease, so that the 
branching does not remain macroscopic in the N —> oo limit, which would conflict with 
the experimental results). It can be seen that the dependence of 8y on A is more or less 
linear. This may be related to the linear increase of the size of the average microbranch as 
a function of the driving displacement (for example, see Fig. [8](b))below). 
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FIG. 3. (Color online) (a) The scaled width of the microbranching pattern &y/w for (a) a honeycomb 


perturbed lattice, (b) the CRN model and (c) a hexagonal perturbed lattice, for different lattice 


sizes. 


IV. RESULTS 


A. The crack’s tip shape 

As mentioned in the introduction, the new experiments using gels yield clear snapshots 
of the crack tip pnun]. Larger scale simulation enable us to compare the crack tip shape 
to the experimental snapshots in both steady-state cracks and near the origin of instability. 

First, we present the very good agreement between experiment and simulations of the 
crack’s tip shape in steady-state cracks in the finite-width strip geometry in Fig [4j In the 
upper picture, when the crack length is small comparing to the sample’s width, the crack 
has the (well-known) parabolic shape; no blunting can be seen and LEFM works perfectly 
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FIG. 4. (Color online) Several snapshots of the crack’s tip in steady-state mode-I fracture in 
experiments on gels. In the upper picture, when the crack length is small comparing to the sample’s 
width, the crack has the (well-known) parabolic shape. As the crack length grows, the finite size 
of the strip effects the crack’s tip shape, generating a “tadpole-like” shape. The lattice simulations 
for the crack tip shape are added to the snapshots in the dashed curves. The experimental pictures 
are taken from f3]. 

(except for nonlinear effects in the extreme crack tip zone [7]). As the crack length grows, 
the finite width of the strip affects the crack’s tip shape, generating a “tadpole-like” shape. 
We set a careful lattice simulations in finite width strip using different values of A/Aq- 
Since the simulations have only a finite number of atoms, the results are scaled to the real 
size of the experimental sample. We can see in Fig. [4] the excellent match of the crack’s 
tip shape between pure lattice simulations (described in detail below) and the experiments. 
This “tadpole” shape is generic for finite-size (strip) lattice simulations, both honeycomb 
or hexagonal, with any amount of viscosity. Since the experimental crack is not exactly 
symmetric between the upper and the lower size of the crack, we used several values of 
A/A g to get an optimal match between experiment and simulation (the small A/ Aq shape 
is somewhat better on the upper side and the large A/ Aq shape is somewhat better in the 
lower side). 
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FIG. 5. (Color online) Several snapshots of the crack tip near the origin of the microbranching 
instability, on the right, there is a series of experimental snapshots of the crack tip before, during, 
and right after a new microbranch is born. On the left, there is a series of snapshots of our 
lattice simulations. Qualitatively, there is a very close resemblance between the two series. The 
experimental pictures are taken from [3]. 

Moreover, in [3] there are snapshots of the crack tip shape near the origin of the micro¬ 
branching instability (on the right in Fig. [5]). In order to reproduce this crack tip shape 
(in addition to using the large scale simulation), we set e = 2ao (only for this part in the 
paper), to magnify the crack tip radius compared to the lattice scale do- This enables us to 
meaningfully compare the simulation results to experiments (this value of e is a bit extreme, 
but qualitatively, the same physical effect can be seen using smaller value of e). The lattice 
simulations are presented on the left in Fig. [5]). 

We can see that there is a close resemblance between the two series’s of snapshots. At 
first, the crack travels in the midline of the sample, yielding the parabolic crack’s tip shape of 
LEFM. When a new microbranch is born, the crack tip bifurcates to two, while very rapidly 
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one “wins” and continues to propagate, and the other “dies”. The branch that continues 
to propagate has a small deformed shape comparing to the original shape before it continue 
to propagate (the crack shape here significantly differs from the LEFM prediction). All 
of these features are reproduced in the lattice simulations. It is important to mention that 
both experiments and simulations reveal that the origin of microbranching in mode-I fracture 
always lies in the immediate vicinity of the crack tip itself, rather than far behind the main 
crack. 


B. Length of microbranches 

In this section we present the quantitative results for various features of the fracture 
patterns, especially the microbranches. First, we present the v(A) curve, the total amount 
of broken bonds in the microbranches, measured for a crack that reached the end of the 
strip, as a function of the crack velocity and the RMS of dR(t)/dt , the rate of increase of 
the derivative of the electrical resistance with respect to time as a function of the crack 
velocity. The results are normalized to to cr, the Rayleigh wave speed, which is calculated 
for different values of kg in Appendix [Cj The results for the perturbed honeycomb lattice is 
presented in Fig. [6j and for the CRN in Fig. [7} 

We can clearly see that in all three models the slope of the r;(A) appears to saturate in 
the high velocity regime for the larger system sizes. This saturation is known from previous 
lattice studies [20]. The shape of the curves are similar to the experimental v(A) curves 
of real amorphous materials [I]. A close look at the curves of the total amount of broken 
bonds in the microbranches (which is our proxy for the average length of a microbranch 
as measured experimentally, which we use to reduce the statistical noise) as a function 
of the crack velocity, using both the honeycomb perturbed lattice and the CRN, reveals 
quantitatively what we have seen qualitatively using small system sizes. Due to the noisiness 
of using finite size lattices, rather than two different regimes, with a single steady state crack 
at small velocities, and a sharp increase in the length of a microbranch in the large velocity 
regime, we get a smooth exponential behavior. However, the exponent of the curves increases 
robustly with the lattice size, sharpening the difference between the steady state regime and 
the the microbranching (large velocity) regime. In addition we can see that in general, 
the CRN results act like a very perturbed honeycomb lattice, i.e., in any given A/Aq, the 
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FIG. 6. (Color online) (a) The v (A) curve of a perturbed honeycomb lattice using different lattice 
sizes, (b) Total number of broken bonds in the microbranches as a function of the crack velocity, 
(c) The RMS of the derivative of electrical resistance with respect to the time as a function of the 
crack velocity. 
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FIG. 7. (Color online) (a) The w(A) curve of a CRN lattice using different lattice sizes, (b) Total 


amount of broken bonds in microbranches, as a function of the crack velocity. 
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FIG. 8. (Color online) (a) The v (A) curve of a perturbed hexagonal lattice using different lattice 


sizes, (b) Total size of broken bonds in microbranches, as a function of the crack velocity. 


number of broken bonds is larger in the CRN than in the perturbed honeycomb lattice, 
similar to what we saw in Ref. [T3]. In addition, using the honeycomb perturbed lattice 
we can see the increase of the RMS amplitude of the derivative of the electrical resistance 
with respect to the time as a function of the crack velocity. At high velocities, the RMS 
amplitude is approximately three times greater than the amplitude in low crack velocities, in 
agreement with the experiments regarding the RMS of the crack velocity (Fig. 11(c) in [5j). 

Looking at the results for the perturbed hexagonal lattice, we get one of the most impor¬ 
tant results of this paper. In Fig. [8j(b) , we present the normalized number of microbranches 
broken bonds (per number of bonds in the entire system) as a function of the crack velocity. 
The different curves (using various lattice size) are similar to each other, however, enlarg¬ 
ing the lattice size, the transition become sharper with /. Using / = 8, we can see two 
different separate regimes (one of steady-state cracks, and one for microbranching), that 
looks very much like the transition of the average microbranch length in the experiments 
(Fig. 11(a) in [5]). This result verifies the main assumption of the lattice models. The 
physical phenomena of microbranching is qualitatively described by lattice models and sim¬ 
ulations, when enlarging the system size, the results become more quantitatively similar to 
the (macroscopic) experimental results. 

Nevertheless, quantitatively, the critical velocity that was found in the experiments was 
less than the critical velocities that were seen in our simulations. However, in previous 
studies it was shown that in the lattice models, we can control the critical velocity using 
different values of a pot [20], I2T] . Here we check that this is still valid when using a 3-body 
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FIG. 9. (Color online) (a) The total size of microbranches as a function of the driving displacement 
A/A q using different values of a po t along with the piecewise linear model (a po t —> oo) using 
perturbed honeycomb lattice with / = 1. (b) The microbranching pattern in a perturbed hexagonal 
lattice using larger scale lattices (/ = 8) for A/Ag = 1.6 and r] = 0.25 for a pot —> oo in the upper 
curve, for a pot = 5 in the lower curve. 

force-law and a perturbed lattice. In Fig. |9](a), for the case / = 1, we present the total 
amount of microbranches using the perturbed honeycomb lattice. We can see that first, 
a pot = 1000 reproduces the piecewise linear results, as expected. Second, smaller values of 
ctpot yields a much larger quantity of broken bonds for a given A, indicating a much lower 
critical velocity (because of the noise in this lattice we cannot realize the zero microbranch 
regime). In addition we check whether using larger-scale lattices (/ = 8) yields the same 
fracture patterns. In Fig. |9](b) we can see that fracture patterns remain similar using the 
smaller value of a pot (here we used the perturbed hexagonal lattice), though for a given A 
the microbranches are larger, as expected. 


C. Microbranching statistics 

Having larger systems enables, for the first time, the generation of enough statistics to 
examine important quantitative features of the microbranches that have been measured 
experimentally Since in this study, we obviously are restricted to 2D features only, 

we focused here on the branching angle of the microbranches, and on the power-law shape 
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FIG. 10. (Color online) The distribution of the branching angle of the microbranches arise nearby 
the main crack using (a) a perturbed honeycomb lattice, (b) the CRN model and (c) a perturbed 
hexagonal lattice. 


of the microbranches. 

The experimental studies on PMMA finds a narrow distribution of the branching angle 
between 20° ^ 6 ^ 40°, with an average angle of 30° 0. We note that a previous study 
using a random perturbed Born-Maxwell model [40] yielded the wrong branching angle, 


namely 15°. In Fig. 10 we present the branching angle distribution of all the microbranches 
generated using all the values of A /Ac in the different models that were used in this study. 
We can see that in all the models studied, the average branching angle is near 30°, very 
much like the experimental branching angle. The variance is different using the different 
models, where the variance of the CRN lattice is the narrowest, and thus most similar to the 
experiments. We also note that this result corresponds nicely with the LEFM microbranch 
analysis of mi which yields a 27° branching angle. 

One of the most striking features that was discovered in the experiments was a power-law 
shape of the branches, y = ax a (x and y are the spatial coordinates of the microbranch such 
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that x = y = 0 is the location of the “root” of the microbranch close to the main crack), 
with a universal power, a ~ 0.7 for several different materials Dig. On the one hand, our 
previous studies using CRN or a perturbed honeycomb lattice yielded only straight (a ~ 1) 
microbranches B3U25]. On the other hand, using a perturbed hexagonal lattice ra we got 
branches that showed a power-law shape with a ~ 0.5, though the results were noisy due to 
the relatively small number of broken bonds in each microbranch. We note that the LEFM 
analysis of [4IJ yields a = 2/3, which is close to the experimental result, but this analysis 
assumes symmetric branching, when our branching (and the experimental results) is by no 
means symmetric. We also note that the branches seen in the elastic beam model [32j or the 
Born-Maxwell model m look very noisy, without a clear power-law shape. In this study, 
using larger lattices we can check the shape of the microbranches based on relatively large 


microbranches 100 broken bonds in each microbranch). In Figs. 11 and [l2[a) we can see 
the power distribution for the different kinds of lattices. 

We can clearly see that the maxima of all the distributions are around a ~ 0.85 — 0.9, 
which is indeed less than 1 (straight lines). Looking closely at the larger lattices (/ = 9) in 
Figs. [T| and [2] we can see the nonlinear power law shape of the large branches, especially for 


the perturbed hexagonal lattice. In Fig. 12 b) we can see a scatter-plot of a as a function of 
the branch size for the perturbed hexagonal lattice case. It is clear that a converges on the 
value a ~ 0.7 for the largest microbranches, which is close to the experimental value. The 


data for all three models is shown in Fig. [13j The other models do not show the decrease 
in a for the longest branches that is apparent for the hexagonal case, and perhaps data for 
longer microbranches are needed in the other cases, but in any case the asymptotic value of 
the exponent in all cases is less than unity. 


V. DISCUSSION 

We have used relatively large lattices fracture simulations using GPU parallel computing 
(with 0 (5 • 10 6 ) particles) for studying mode-I fracture. We find that the basic results 
from small lattices (0(5 ■ 10 4 ) particles) are confirmed using the larger size systems. The 
fracture patterns look more physical with a larger system, due to the large number of 
broken bonds in each microbranch. The width of the microbranch region relative to the 
system width decreases, a necessary condition if these models are to be taken seriously. 
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FIG. 11. (Color online) The distribution of the power of the power-law shaping of the microbranches 
arise nearby the main crack using (a) a perturbed honeycomb lattice and (b) the CRN model. 
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FIG. 12. (Color online) (a) The distribution of a, the power-law exponent (y = ax a ) of the 
nricrobranches arising near the main crack for a perturbed hexagonal lattice, (b), A scatter-plot 
of the power a of the power-law exponent of the microbranches. We can clearly see that for the 
long branches, a converges on the value a ~ 0.7, which is close to the experimental value. 



The basic properties of the nricrobranches, like the total length of the nricrobranches and 
the oscillations of the derivative of the electrical resistance with respect to the time based 
velocity measurements, which were extremely noisy when using small lattices, look more 
smooth and realistic in the larger lattices. In particular there is now a clearer transition 
between the regime of steady state cracks, which there is a negligible amount of broken bonds 
beside the main crack, and the regime of instability, where the amount of broken bonds in 
the microbranches increases dramatically. The sharp transition is particularly clear in the 
hexagonal perturbed lattice (Fig. |8](b)). 

In addition, important features of the nricrobranches that have been found and studied 


21 


































































1.5 



FIG. 13. (Color online) The distribution of a , the power-law exponent (y = ax a of the micro¬ 
branches arising near the main crack for the perturbed hexagonal, perturbed honeycomb and CRN 
models, as a function of the micro branch length. 

experimentally, are recovered in our large system simulations. The correct branching angle 
is found, and in the CRN lattice the correct variance is also obtained. The universal power- 
law shape that was found in different experimental studies [5J is recovered here, and for the 
larger branches, we get the correct power in the hexagonal perturbed lattice. 

In future work, we plan to exploit the power of GPU parallel computing to run 3D simu¬ 
lations using 0(5 ■ 10 6 ) particles, with the goal of studying the different aspects concerning 
the 3D nature of the microbranches. We intend to check the similarities and the differences 
between 2D and 3D Mode-I fracture simulations, and to find the regime when the 2D model 
is sufficient, and on the other hand, the regime where 3D simulations are crucial. 


Appendix A: GPU-accelerated C Code 

In this appendix we discuss implementation of our codes using CUDA and the run-times 
of the CUDA runs using NVIDIA’s Tesla C2050/C2070 GPU for the different modules that 
were used in this study, and the acceleration ratios between a single CPU run and the GPU 
run. We re-wrote our C codes and replaced the main time-demanding functions by CUDA 
kernels Ed- CUDA is a computer code language that was developed by NVIDIA for using 
its CPU’s. We note that in principle, one could use OpenCL, a general language for any 
GPU device, though CUDA is optimized for NVIDIA’s graphic cards. 
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For using GPU computing optimally, the code should be written as a function of the 
particular hardware. Naively, the computational tasks are divided into several blocks, each 
block containing a specific number of threads, when the threads in the same block execute 
simultaneously. In this work we used 512 threads per block which is the maximal number of 
threads per block available on our graphic card. The parameters of the physical problem are 
loaded to the global memory of the GPU and each thread computes one computational task, 
such as the force of one spring that acts on a pair of atoms. Such a naive choice however only 
produces up to a 10 times acceleration, since there are a lot of calls to the global memory 
(which is relatively slow) for each atom calculating the force, especially in the 3-body force- 
law. Thus we extensively employed the option of using shared memory (which is as fast as 
cache memory, shared for all the threads in the same block and limited to 64K) to further 
accelerate the simulations, especially in the molecular-dynamics module [371 38] . The global 
memory of Tesla C2050/C2070 GPU is about 2.5 Gbytes, which is the limiting factor of the 
number of atoms in the simulation. In this work we used double precision accuracy for all 
of our calculations (so using float precision will increase the number of atoms by a factor of 
2, which is not significant; the main issue of this study is the effect of increasing the number 
of atoms by two orders of magnitude). 

The underlying plan of the main, molecular dynamics, module using shared memory is 
to split the (potentially random, and thus general) grid of atoms to several physical zones 
(very much like Open-MPI as opposed to OpenMP parallelization) with lists that connect 
between the zones applied as a boundary condition for each physical zone; for the central 
elastic force law (and also for the viscoelastic force law) we sort the bonds and for the 
3-body force law we sort the atoms. Each zone (“block” in the CUDA-jargon) loads the 
locations and the velocities to a fast shared memory, and each “thread” calculates the force 
of a certain bond (for the central force law), or atom (for the 3-body force law). Thus, 
instead of several calls to the global memory for each atom’s location, the calls are for the 
fast shared memory, and so good efficiency is then achieved. Since the threads in the same 
block execute simultaneously, we have to use CUDA Atomic commands to sum correctly the 
contribution of the forces acts from the neighboring atoms. Then, the new velocities and 
locations are calculated by simple CUDA kernels. 

The electrical resistance module is basically solving a nonlinear Laplace equation. We 
used the methodology that was introduced in [39] for calculating the electrical resistance (a 
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FIG. 14. (Color online) The derivative of the electrical resistance with respect to the time using 
a perturbed honeycomb lattice with a /a g = 3.4 using the different lattice sizes. The plots are 
normalized in the x-axis to the / = 1 size (i.e. the / = 9 results are divided by 9 etc.) 


constant grid of bonds with a = 1), while a cracked bond in the molecular-dynamics module 
determines the “cracked” (er = 0) bonds in the constant grid of the electrical resistance. 
Thus, we used the same well-known methodology of using the CUDA kernels for solving 2D 
diffusion equations [3B|, including the use of shared memory. We implemented both Jacobi 
and red-black gauss-Seidel methods of solution, but no significant difference (in terms of the 
number of iterations for convergence) was found between the methods. In Fig. 14 we can 


see an example of the derivative of the electrical resistance with respect to the time using a 
perturbed honeycomb lattice for a specific a /a g . We can see that using different size lattices 
(using the new GPU code in the larger lattices and the CPU code in the smaller lattice), 
we can see that the shape of the curves of electrical resistance look very much alike the 
experimental RMS amplitude of the crack velocity [5J [6] (that is measured via the electrical 
resistance ESI)- 

The parallel Monte-Carlo algorithm module for generating the CRN required extra care. 
Since each possible switch of bonds should be considered energetically independently, each 
switch should be separate from all the other simultaneous possible switches (to avoid over¬ 
lapping of the switches and their neighboring zone). We mention that we use the parallel 
THRUST library (in CUDA) for sorting efficiently the nearest bonds for each bond, every 
given number of cycles. 


In Figs. 15 16](a) we can see the typical run times for 1000 cycles (in seconds) for the 
molecular dynamics module and the Laplace solver module as a function of the system size. 
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FIG. 15. (Color online) (a) Simulation times (in [sec]) of 1000 cycles of the molecular-dynamics 
Euler scheme as a function of the system size (both using honeycomb and hexagonal lattices) 
using unoptimized C code with CPU and with CUDA using GPU. (b) The acceleration run times 
between GPU and CPU as a function of the system size for honeycomb and hexagonal lattices). 


We can clearly see the benefit of using GPU computing, while the main benefit is the possi¬ 
bility to run systems with large number of atoms, which with a single CPU, was prohibitively 
time-consuming. The run-time using (9(10 4 ) atoms with a single CPU is similar to the run¬ 


time using 0( 10 6 ) atoms with a single GPU. In Figs. 15 16(b) we see the acceleration ratios 
between a single CPU and a single GPU (of course, only for small systems, when a CPU 
run is available). We can see the significant acceleration, 40 times faster for the honeycomb 
lattice and over 50 times faster for a hexagonal lattice due to more demanding 3-body force 
law). In the Laplace solver the speedup is a little bit lower and stands at approximately 25 
times faster in GPU versus a single CPU. 


In Fig. [17] we can see the run-times using a single CPU and the GPU. We can see here 
that the acceleration ratio here is lower then in the previous modules (about « 5 - 10), 
but still, in larger lattices (of 0( 10 6 ) atoms), the benefit is clear. We mention that since 
the programing using CUDA is much more demanding from the programmer, especially 
regarding the memory management, while re-programing the code, we improved our old- 
CPU code (that was in use in [25]), significantly; the acceleration ratio of the GPU code to 
the old CPU code is ^60. 
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FIG. 16. (Color online) (a) Simulation times (in [sec]) for 1000 iterations of a Jacobi-method 
Laplace solver (for the electrical resistance) as a function of the system size using unoptimized C 
code with CPU and with CUDA using GPU (the times are similar also in Red-Black Gauss Seidel 
method), (b) The acceleration run times between GPU and CPU as a function of the system size. 




FIG. 17. (Color online) (a) Simulation times (in [sec]) for a WWW Monte-Carlo scheme for 
producing a CRN (simulation stopped while the energy reaches half of the initial random energy) 
as a function of the system size using unoptimized C code with CPU, with CUDA using GPU and 
with the “old-CPU” code used in [25] . (b) The acceleration run times between GPU and CPU and 
the “old-CPU” codes as a function of the system size. 


Appendix B: CRN Monte-Carlo parallel CUDA algorithm 


In previous work ra we have shown that the CRN shares similar features with real 
amorphous matter, like amorphous silicon [32] . In this appendix we check explicitly the 


quality of the parallel GPU algorithm, generating a CRN. In Fig. 18 we can see the radial 
and angular distributions of the bonds in the CRN using different lattice size, and in Fig. [19 
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FIG. 18. (Color online) (a) The radial distance distribution of the bonds in CRN using different 
lattice size, (b) The angular distribution of the bonds in the CRN using different lattice size. 



FIG. 19. (Color online) The radial distribution function (RDF or g(r)) of the the CRN using 
different size of lattices. 


we can see the radial distribution function g{r ) using different lattice size. / = 1 (of 0( 10 4 ) 
atoms) is the data using the old CPU code that was in use in [25J, when larger lattices was 
produced via the new GPU algorithm. 

We can see that the radial and angular distributions and the g(r) curves are similar under 
the scaling of the lattice size (the number of the bonds or angles). This proves the validity 
of the parallel GPU algorithm, comparing the old CPU algorithm (that was verified before 
against experiments). Moreover, the RDF in the larger lattice size is smoother due to better 
statistics. 
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Appendix C: The Rayleigh surface wave speed with kg ^ 0 lattices 


Since the models in this paper use a 3-body potential law (aside by the central two-body 
force law) with k g % 0, we need to recalculate the Rayleigh wave speed cr , which is the 
terminal velocity for mode-I fracture pj for different kg/k r . The most convenient way to 
calculate the Rayleigh wave speed is to calculate first the longitude (primary) q and the 
transverse (secondary) q wave speeds and then, to calculate the Rayleigh wave speed via 
the well-known formula [ 45] : 


1 - 




(Cl) 


For kg = 0 the expressions for q and q, and thus also for cr, can be derived analytically for 
a 2D hexagonal lattice with a as the lattice scale (for k r — rn — 1; the wave velocities scales 
as yjk r /m)\ 
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and for a 2D honeycomb lattice: 
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We calculate q and q via measuring the wave velocities by initiating longitude and 
transverse small deformations in the end of the samples in the different lattices that we use 


in this study and then find cr via Eq. Cl The results are shown in Fig. |20[a) for the 


hexagonal lattice and in Fig. 20 b) for the honeycomb and CRN lattices. 

We can see that for both lattices, the numerical value for the wave velocities using kg = 0 


reproduce the analytical values, Eqs. |C2| and [C3| respectively. In the hexagonal lattice using 
kg/k r = 10 (which was the value used in this study) the Rayleigh wave speed increases by 
~ 65% relative to the kg — 0 value. I 11 the honeycomb lattice, where we use kg/k r = 1 at 
most, the Rayleigh wave speed increases by ~ 12% relative to the kg — 0 value. We note the 


CRN speeds (in the ellipse in Fig. 20 b)) are a little bit slower than the pure honeycomb 
lattice. I 11 addition we note that the random noise of the perturbed lattice changes the wave 
speeds less than 1%, and that the wave speeds are not affected at all by a po t- 











(a) 



FIG. 20. (Color online) (a) The longitude and the transverse sound waves speeds along with 


the resulting calculated Rayleigh surface wave speed using Eq. Cl for the hexagonal lattice as a 


function of kg/k r . (b) The same in the honeycomb lattice. The sounds velocities for the CRN that 
was used in this paper (ko/k r = 1) are presented In the ellipse. 
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